#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__1 = 1;


/* ------------------------------------------------------------------ */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dortqr_
   (integer *nz, integer *n, integer *nblock, doublereal *z__, doublereal *b)
{
  /* System generated locals */
  integer z_dim1, z_offset, b_dim1, b_offset, i__1, i__2;
  doublereal d__1;

  /* Builtin functions */
  //double d_sign(doublereal *, doublereal *);

  /* Local variables */
  /* extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
  integer *);*/
  static doublereal temp;
  /* extern doublereal dnrm2_(integer *, doublereal *, integer *);*/
  static integer i__, j, k, m;
  //extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
  // integer *);
  static doublereal sigma;
  //extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
  // integer *, doublereal *, integer *);
  static integer length;
  static doublereal tau;



  /* THIS SUBROUTINE COMPUTES THE QR FACTORIZATION OF THE N X NBLOCK */
  /* MATRIX Z.  Q IS FORMED IN PLACE AND RETURNED IN Z.  R IS */
  /* RETURNED IN B. */


  /* THIS SECTION REDUCES Z TO TRIANGULAR FORM. */

  /* Parameter adjustments */
  z_dim1 =  *nz;
  z_offset = 1+z_dim1 * 1;
  z__ -= z_offset;
  b_dim1 =  *nblock;
  b_offset = 1+b_dim1 * 1;
  b -= b_offset;

  /* Function Body */
  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {

    /* THIS FORMS THE ITH REFLECTION. */

    length =  *n - i__ + 1;
    d__1 = dnrm2_(&length, &z__[i__ + i__ * z_dim1], &c__1);
    sigma = d_sign(&d__1, &z__[i__ + i__ * z_dim1]);
    b[i__ + i__ * b_dim1] =  - sigma;
    z__[i__ + i__ * z_dim1] += sigma;
    tau = sigma * z__[i__ + i__ * z_dim1];
    if (i__ ==  *nblock)
    {
      goto L30;
    }
    j = i__ + 1;

    /* THIS APPLIES THE ROTATION TO THE REST OF THE COLUMNS. */

    i__2 =  *nblock;
    for (k = j; k <= i__2; ++k)
    {
      if (tau == 0.)
      {
        goto L10;
      }
      temp =  - ddot_(&length, &z__[i__ + i__ * z_dim1], &c__1, &z__[i__ + k *
        z_dim1], &c__1) / tau;
      daxpy_(&length, &temp, &z__[i__ + i__ * z_dim1], &c__1, &z__[i__ + k *
        z_dim1], &c__1);
      L10: b[i__ + k * b_dim1] = z__[i__ + k * z_dim1];
      z__[i__ + k * z_dim1] = 0.;
      /* L20: */
    }
    L30: ;
  }

  /* THIS ACCUMULATES THE REFLECTIONS IN REVERSE ORDER. */

  i__1 =  *nblock;
  for (m = 1; m <= i__1; ++m)
  {

    /* THIS RECREATES THE ITH = NBLOCK-M+1)TH REFLECTION. */

    i__ =  *nblock + 1-m;
    sigma =  - b[i__ + i__ * b_dim1];
    tau = z__[i__ + i__ * z_dim1] *sigma;
    if (tau == 0.)
    {
      goto L60;
    }
    length =  *n -  *nblock + m;
    if (i__ ==  *nblock)
    {
      goto L50;
    }
    j = i__ + 1;

    /* THIS APPLIES IT TO THE LATER COLUMNS. */

    i__2 =  *nblock;
    for (k = j; k <= i__2; ++k)
    {
      temp =  - ddot_(&length, &z__[i__ + i__ * z_dim1], &c__1, &z__[i__ + k *
        z_dim1], &c__1) / tau;
      daxpy_(&length, &temp, &z__[i__ + i__ * z_dim1], &c__1, &z__[i__ + k *
        z_dim1], &c__1);
      /* L40: */
    }
    L50: d__1 =  - 1. / sigma;
    dscal_(&length, &d__1, &z__[i__ + i__ * z_dim1], &c__1);
    L60: z__[i__ + i__ * z_dim1] += 1.;
    /* L70: */
  }
  return 0;
} /* dortqr_ */
